%% 脚本说明
%{
- --------------------------------------------------
* 文件名:HTL
* 函数名:None
* 功 能:通过傅里叶变换，计算HTL的介电常数和折射率等
* 变量说明:
** Resut1:直接导入两个信号的比值
* 状态:程序已经写完但是没有得到想要的结果，目前直接采用别人的结果
* 注意事项:
- --------------------------------------------------
%}
%% Data Preparation
d = 823e-6; % thickness
c = 3e8; % speed of light

%%
% input the ratio of the sample and reference signals
Result = Resut1;
% w = (0:(2048)/2)/2048*1e12; % w is the freq
w = w1 * 1e12;
Phi = -angle(Result);
% n = 1;

unphasewrap = Phi;

for i = 2:length(unphasewrap)
    diff = unphasewrap(i) - unphasewrap(i - 1);

    if diff >= pi
        unphasewrap(i:end) = unphasewrap(i:end) - 2 * pi;
        %         n = n + 1;
    elseif diff <= -pi
        unphasewrap(i:end) = unphasewrap(i:end) + 2 * pi;
    end

end

% Phi = unwrap(Phi);
Phi = unphasewrap;

ns = Phi ./ w * c / d + 1;
ks = log(4 * ns ./ abs(Result) ./ (ns + 1).^2) ./ w * c / d;
ks = real(ks);
alphas = 2 * ks .* w / c;
epsilon1 = ns.^2 - (c * alphas / 2 ./ w).^2;
epsilon2 = 2 .* ns .* ks;
epsilon = epsilon1 - 1j * epsilon2;

figure
subplot(511);
plot(w, ns)
legend('n'); xlim([0.2e12 4e12]);
subplot(512);
plot(w, ks);
legend('A'); xlim([0.2e12 4e12]);
subplot(513);
plot(w, (epsilon1));
% ylim([0 3]);
xlim([0.2e12 4e12]);
legend('\epsilon1');
subplot(514);
% plot(w, angle(epsilon));
plot(w, tan(epsilon2 ./ epsilon1));
legend('tan');
% ylim([-3 3]);
xlim([0.2e12 4e12]);
subplot(515)
plot(w, Phi)
legend('Phi')
xlim([0.2e12 4e12]);
